
% Driver script for the function computererrbounds_network()
% 
% Anand Sengupta and P. Ajith, 8 Aug 2009
% 
% $Id:$

 clear all;

 setconstants;

 MVec		= logspace(log10(2), log10(40), 15);	% in Solar mass
 detNameVec	= {'I1'};
 noisePSDVec	= {'AdvLIGO'};
 lambdaVec	= {'M', 'eta', 'tc'};
 %lambdaVec	= {'a1', 'a2', 'a3', 'a4', 'tc'};

 Params.isAngAvged	= 0;
 Params.tc		= 128;			% Coalescence time			
 Params.fLower		= 30.0;			% Lower cut-off frequency
 Params.eta		= 0.25;			% symmetric mass ratio
 Params.Alpha		= 1.354957275964931*rand(1,1);	% Overhead on Bangalore
 Params.Delta		= 0.226311026341931*rand(1,1);	% Overhead on Bangalore
 Params.psi		= 0.7824675;			% Polarization angle
 Params.startPhase	= 0.6757834;
 Params.iota		= 0;			% inclination angle
 Params.r		= 100*1e6*PARSEC_SEC;	% Mpc
 Params.gmstRad		= 0;
 Params.A		= [];
 Params.B		= [];
 Params.P		= [];
 Params.Q		= [];

 %--- Calculate the derived parameters
 c	= cos(Params.iota);
 Phi0	= Params.startPhase;
 psi	= Params.psi;

 Params.a1 = ( -(1+c^2)*cos(Phi0)*sin(2*psi)...
       - 2*c*cos(2*psi)*sin(Phi0) ) / Params.r;

 Params.a2 = ( 2*c*cos(Phi0)*sin(2*psi)...
       + (1+c^2)*cos(2*psi)*sin(Phi0) ) / Params.r;
   
 Params.a3 = ( -(1+c^2)*cos(Phi0)*cos(2*psi)...
       + 2*c*sin(2*psi)*sin(Phi0) ) / Params.r;

 Params.a4 = ( -2*c*cos(Phi0)*cos(2*psi)...
       + (1+c^2)*sin(2*psi)*sin(Phi0) ) / Params.r;

 %--- Clear temporary variables
 clear c Phi0 psi

 % compute the error bounds for each pair of masses
 for iM=1:length(MVec)

    Params.totalMass = MVec(iM);

    [sigma(iM,:), covMat, snr(iM), snrVec, condNum(iM), condNumVec, Params ] = ...
    	computeerrbounds_network(Params, detNameVec, noisePSDVec, lambdaVec);
    
 end


 %--- Plot the RMS errors against the total mass
 for iParam = 1:length(sigma(1,:))

    if ( sum(strcmp(lambdaVec{iParam}, {'M', 'eta', 'tc'})) )
        x = ones(size(MVec));
    else
	x = getfield(Params, lambdaVec{iParam})*ones(size(MVec)) 
    end
       
    figure(iParam); clf;

    %-- Percentage error
    y = abs( sigma(:,iParam)*100 ./ x(:) );
    
    semilogx(MVec, y,...
		'ro-', 'MarkerFaceColor', 'b');

    if (max(sigma(:,iParam))/min(sigma(:,iParam)) >= 6)
	set(gca, 'YScale', 'log'); 
    end	 

    grid on; axis tight;
    set(gca, 'FontSize', 12);
    set(gca, 'Xlim', [0.9 1.4*max(MVec)]);
    xlabel('M / M_{sun}', 'FontSize', 14);
    ylabel(sprintf('\\Delta %s', lambdaVec{iParam}), 'FontSize', 14);

    clear x;

 end

